The 4th Dimension of Movement

Analysis (empirical tracks)

Authors

Johannes Signer

Cédric Scherer

Stephanie Kramer-Schadt

Published

February 26, 2025

Preparation

library(readr)
library(dplyr)
library(tidyr)
library(forcats)
library(stringr)
library(ggplot2)
library(ggtext)
library(lubridate)
library(purrr)
library(sf)
library(terra)
library(rnaturalearth)
library(rcartocolor)
library(colorspace)
library(patchwork)
library(amt)  
library(corrr)
library(here)

## plot theme
theme_set(d6::theme_d6(base_size = 18))

theme_update(
  panel.grid.major.x = element_blank(),
  panel.grid.major.y = element_line(color = "#fefefe", linewidth = .5),
  panel.grid.minor = element_blank(),
  panel.background = element_rect(fill = "grey95", color = "transparent"),
  panel.border = element_rect(fill = "transparent", color = "transparent"),
  strip.text = element_markdown(color = "black", hjust = .5, face = "bold", size = rel(1.3)),
  legend.text = element_text(family = "PT Mono"),
  legend.position = "top",
  legend.box = "vertical",
  legend.spacing = unit(.2, "lines"),
  axis.text = element_text(family = "PT Mono"),
  axis.text.x = element_text(margin = margin(t = 5)),
  axis.text.y = element_text(margin = margin(r = 5)),
  axis.ticks.length = unit(0, "lines"),
  panel.spacing.x = unit(.7, "lines"),
  panel.spacing.y = unit(.9, "lines"),
  plot.margin = margin(rep(.5, 4))
)

## path to processed data
dir_proc <- here("output", "empirical", "data_proc")

## path to analysis results
dir_est <- here("output", "empirical", "results")

## path + folder for movement plots
dir_move <- here("plots", "emp_tracks")
if (!dir.exists(dir_move)) dir.create(dir_move, recursive = TRUE)

## path + folder for plots
dir_plots <- here("plots")
if (!dir.exists(dir_plots)) dir.create(dir_plots, recursive = TRUE)

Data Import

Movement Data

df_boars <- read_rds(paste0(dir_proc, "/boars_resampled_3035.rds"))

sf_boars <- 
  df_boars %>% 
  filter(!is.na(x), !is.na(y)) %>% 
  st_as_sf(
    coords = c("x", "y"), 
    crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
  ) %>% 
  st_make_valid() %>% 
  st_transform(crs = 3035) %>% 
  mutate(
    x = st_coordinates(.)[,1],
    y = st_coordinates(.)[,2],
    t = ymd_hms(time),
    id = str_replace(id, "_", " "),
    res_label = paste(res, "minutes"),
    res_label = fct_reorder(res_label, res)
  )

## color palette for boars
n <- length(unique(df_boars$id))
pal <- carto_pal(name = "Bold", n = n+1)[1:n]
names(pal) <- unique(df_boars$id)

Landscape Variables

## turn into data frame for plotting
ras_env <- rast(here("output", "empirical", "data_proc", "env_vars_100m.tif"))

df_env <- 
  ras_env %>%
  terra::as.data.frame(xy = TRUE) %>% 
  as_tibble() %>% 
  mutate(hab = factor(hab, levels = c(0, 1), labels = c(" non-forest   ", " forest")))

bln <- sf::st_transform(d6berlin::sf_berlin, crs = st_crs(ras_env))

Figure 2

Movement tracks of the wild boars

sf_boars_move_12h <- 
  sf_boars |> 
  filter(res == 720) |>
  group_by(id) |>  
  mutate(t_id = row_number()) |>
  arrange(time)
  
sf_boars_centr_12h <-
  sf_boars_move_12h %>% 
  group_by(id) %>%
  summarize(geometry = st_union(geometry)) %>%
  st_centroid() %>%
  mutate(x = st_coordinates(.)[,1], y = st_coordinates(.)[,2])

set.seed(1)

## Plot
p_1 <- 
  ggplot(df_env) + 
    geom_tile(aes(x, y, fill = imp)) + 
    geom_sf(data = bln, fill = NA, linewidth = 2.7, color = "white") +
    geom_sf(data = bln, fill = NA, linewidth = .9, color = "#415c27") +
    ggforce::geom_mark_hull(
      data = sf_boars_centr_12h,
      aes(x, y, group = id, label = paste0("WB", id)), 
      color = "transparent", con.colour = "grey30", 
      label.colour = "grey20", label.family = "PT Sans Narrow", label.fontsize = 20,
      label.buffer = unit(12, 'mm'), label.margin = margin(2, 2, 1, 2, 'mm'), 
      expand = unit(1, 'mm'), con.cap = unit(.01, 'mm'), con.size = .7
    ) +
    geom_sf(
      data = sf_boars_move_12h, 
      aes(group = id),
      inherit.aes = FALSE,
      color = "black",
      size = 1.3
    ) +
    geom_path(
      data = sf_boars_move_12h, 
      aes(x, y, group = id),
      inherit.aes = FALSE,
      linewidth = .8
    ) +    
    geom_sf(
      data = sf_boars_move_12h, 
      aes(color = t_id, group = id),
      size = .7,
      stroke = .3
    ) +
    scale_color_carto_c(palette = "PinkYl", name = "Movement step (relative)") +
    scale_fill_gradient(low = "grey95", high = "grey40", name = "Imperviousness (scaled)") + 
    coord_sf(expand = FALSE, clip = "off") +
    guides(
      fill = guide_colorbar(
        order = 2,
        barwidth = unit(24, "lines"), barheight = unit(.6, "lines"), 
        title.position = "top", title.hjust = .5
      ),
      color = guide_colorbar(
        order = 1, 
        barwidth = unit(24, "lines"), barheight = unit(.6, "lines"),
        title.position = "top", title.hjust = .5)
    ) +
    labs(x = "Longitude", y = "Latitude") +
    theme(
      plot.title = element_text(hjust = .5), 
      legend.box = "horizontal",
      legend.box.margin = margin(t = -5)
    )

ggsave(paste0(dir_plots, "/2_emp_boars_all_12h.png"), width = 12, height = 8, dpi = 600, bg = "white")

p_1

Single tracks per ID and sampling interval

plot_move_map <- function(boar) {
  
  sf_id <- 
    sf_boars |> 
    filter(id == boar) |> 
    add_count(res) |> 
    mutate(label = paste0(res, " minutes<br><span style='font-weight:400;font-size:18pt;'>(", n, " locations)</span>")) |>
    arrange(res, time) |> 
    mutate(label = fct_inorder(label))
  
  bbox <- 
    sf_id|>
    st_buffer(1000) |> 
    st_bbox()
  
  ## Plot
  g <- 
    ggplot(df_env, aes(x, y)) + 
      geom_tile(aes(fill = imp)) +
      geom_sf(data = sf_id, color = "grey10", size = 1.7, stroke = .3) +
      geom_path(data = sf_id, aes(x, y, group = boar), inherit.aes = FALSE,
                color = "grey10", linewidth = .3) +  
      geom_sf(data = sf_id, aes(color = boar), size = .7) +
      coord_sf(xlim = c(bbox$xmin, bbox$xmax), 
               ylim = c(bbox$ymin, bbox$ymax), expand = FALSE) + 
      facet_wrap(~label, nrow = 1) +
      scale_color_manual(values = pal, guide = "none") +
      scale_fill_gradient(low = "grey95", high = "grey40", guide = "none") +
      labs(x = NULL, y = NULL, title = paste("Wild boar", boar)) +
      theme(
        plot.title = element_text(size = rel(2)),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        panel.spacing.x = unit(2.3, "lines")
      )
  
  ggsave(paste0(dir_move, "/trajectory_boar_", boar, ".png"), 
         width = 25, height = 8.5, dpi = 600, bg = "white")
  
  return(g)
}

if(!file.exists(paste0(dir_move, "/trajectory_boar_9.png"))) {
  walk(unique(sf_boars$id), ~plot_move_map(.x))
}

## example maps
knitr::include_graphics(paste0(dir_move, "/trajectory_boar_5.png"))

knitr::include_graphics(paste0(dir_move, "/trajectory_boar_7.png"))

Figure 6

dat_est_ssf <- read_rds(paste0(dir_est, "/ssf.rds")) |>
  dplyr::select(id, res, fits) |>
  mutate(run_id = 1:nrow(pick(id))) |>
  unnest(cols = fits) |> mutate(method = "issf")

dat_est_rsf <- read_rds(paste0(dir_est, "/rsf.rds")) |>
  dplyr::select(id, res, fits) |>
  mutate(run_id = 1:nrow(pick(id))) |>
  unnest(cols = fits) |> mutate(method = "rsf")

dat_est_wrsf <- bind_rows(
    map(list.files(paste0(dir_est, "/wrsf"), full.names = TRUE), ~
          read_rds(.x) |> mutate(run_id = str_extract(.x, "inst_\\d{1,2}") |> str_extract("\\d{1,2}") |> as.numeric()))
  ) |>
  mutate(term = case_when(term == "habitat" ~ "hab",
                          term == "dist_water" ~ "dtw",
                          term == "imperviousness" ~ "imp"))

dat_est <- 
  bind_rows(
    dat_est_wrsf |> left_join(dat_est_rsf |> dplyr::select(id, res, term, run_id)),
    dat_est_rsf,
    dat_est_ssf
  ) |>
  dplyr::filter(term %in% c("hab", "imp", "dtw")) |>
  mutate(
    id = factor(id, levels = sort(as.numeric(as.character(unique(id))))),
    sim = factor(method, levels = c("issf", "rsf", "wrsf"),
                 labels = c("iSSA", "RSA", "wRSA")),
    c_name = factor(term,
                    levels = c("hab", "imp", "dtw"),
                    labels = c("Forest<br>habitat", "Impervi-<br>nousness",
                               "Distance<br>to water"))
  )

## linear numeric scale with significant outcomes only
dat_est |>
  filter(estimate > -20) |>
  mutate(sig = fct_rev(ifelse(!(low < 0 & high > 0), "significant", "not significant"))) |>
  ggplot(aes(x = res,  y = estimate, group = id)) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0, alpha = .12) +
  geom_hline(aes(yintercept = 0), color = "grey67") +
  geom_line(
    aes(color = id, color = after_scale(desaturate(lighten(color, .33), .5))),
    linewidth = .6
  ) +
  geom_point(
     fill = "white", color = "white", shape = 21, stroke = .6, size = 2.2
  ) +
  geom_point(
    aes(fill = id, alpha = sig), shape = 21, stroke = .6,
    size = 2.2
  ) +
  facet_grid(
    c_name ~ sim,
    scales = "free_y",
    space = "free_y"
  ) +
  coord_cartesian(clip = "off") +
  scale_x_continuous(
    breaks = sort(unique(dat_est$res)),
    labels = sort(unique(dat_est$res)),
    expand = c(.01, .01)
  ) +
  scale_y_continuous(expand = c(.07, .07), breaks = seq(-20, 5, by = 1)) +
  scale_color_manual(
    values = pal, name = NULL, labels = function(x) paste0("WB", x), 
    guide = guide_legend(nrow = 1, order = 1)
  ) +
  scale_fill_manual(
    values = pal, name = NULL, labels = function(x) paste0("WB", x), 
    guide = guide_legend(nrow = 1, order = 1)
  ) +
  scale_alpha_manual(values = c(1, .25), guide = guide_legend(override.aes = list(fill = "grey40"), order = 2)) +
  labs(x = "Sampling interval in minutes", y = "Estimated selection parameter", alpha = NULL) +
  theme(
    axis.text.x = element_text(size = rel(.7), hjust = c(1.05, .25, .25, .4, .5, .75)),
    axis.text.y = element_text(size = rel(.7)),
    axis.title.y = element_text(margin = margin(r = 3)),
    panel.spacing.x = unit(1.4, "lines"),
    legend.position = "bottom",
    legend.text = element_text(family = "PT Sans Narrow"),
    legend.justification = "left",
    legend.box.just = "left",
    #legend.location = "plot",
    legend.key.width = unit(1.4, "lines"),
    legend.margin = margin(1, 1, 1, 1),
    strip.text.x = element_markdown(hjust = 0),
    strip.text.y = element_markdown(angle = 0, hjust = 0, vjust = .5, lineheight = .9,
                                    margin = margin(l = 6, t = -.3)),
    strip.background = element_rect(color = "transparent", fill = "transparent")
  )

ggsave(paste0(dir_plots, "/6_emp_est_env.png"), width = 10.5, height = 11.5, dpi = 600, bg = "white")

Figure S6

Maps of the environmental variables and the study area

Plot Environmental Variables

phab <- 
  ggplot(df_env, aes(x, y, fill = hab, color = after_scale(fill))) + 
    geom_tile(size = .1) +
    scale_fill_manual(
      values = c("#d7b04e", "#005a00"), 
      name = "Forest habitat"
    ) +
    coord_equal(expand = FALSE, clip = "off") +
    guides(fill = guide_legend(title.position = "top", title.hjust = .5)) +
    theme(axis.title = element_blank(), legend.text = element_text(size = 15, family = "PT Sans"))

pimp <- 
  ggplot(df_env, aes(x, y, fill = imp, color = after_scale(fill))) + 
    geom_tile(size = .1) +
    ggsci::scale_fill_material(
      palette = "light-green",
      name = "Imperviousness (scaled)"
    ) +
    coord_equal(expand = FALSE, clip = "off") +
    guides(fill = guide_colorbar(
      barwidth = unit(24, "lines"), barheight = unit(.6, "lines"),
      title.position = "top", title.hjust = .5)
    ) +
    theme(axis.title = element_blank())

pdtw <- 
  ggplot(df_env, aes(x, y, fill = dtw)) + 
    geom_tile(size = .1) +
    scale_fill_carto_c(
      palette = "ag_Sunset",
      direction = -1,
      limits = c(0, NA),
      name = "Distance to water (scaled)"
    ) +
    coord_equal(clip = "off", expand = FALSE) +
    guides(fill = guide_colorbar(
      barwidth = unit(24, "lines"), barheight = unit(.6, "lines"),
      title.position = "top", title.hjust = .5)
    ) +
    theme(axis.title = element_blank())

p_s6a <- 
  (pdtw + labs(x = "", y = "") + 
   pimp + labs(x = "", y = "Latitude") + 
   phab + labs(x = "Longitude", y = "")) * 
  theme(
    plot.margin = margin(5, 15, 5, 0),
    panel.background = element_rect(fill = "grey85", color = "grey85"),
    axis.text = element_text(size = rel(.75)),
    legend.title = element_text(size = rel(1.4), face = "bold")
  ) + 
  plot_layout(ncol = 1)

ggsave(paste0(dir_plots, "/_raw/s6a_maps_emp_env_vars.png"), width = 8.6, height = 18, dpi = 600, bg = "white")

p_s6a

Overview Map

globe <- d6berlin::globe(size_pin = 3, col_earth = "grey96")

sf_europe <- 
  ne_countries(continent = "Europe", scale = 10, returnclass = "sf") |> 
  st_transform(crs = st_crs(sf_boars))

bbox <- 
  sf_boars |> 
  st_bbox() |> 
  st_as_sfc()

germany <-
  ggplot(sf_europe) + 
  geom_sf(aes(fill = sovereignt == "Germany"), color = NA, size = 0) +
  geom_sf(data = bbox, fill = "#415c27", size = 1.2, color = "black") +
  coord_sf(xlim = c(3990000, 4720000), ylim = c(2660000, 3580000), expand = FALSE) +
  scale_fill_manual(values = c("grey92", "#c3dbab"), guide = "none") +
  theme_void()

p_s6b <- globe / plot_spacer() / germany + plot_layout(heights = c(1, .02, 1.6))

ggsave(paste0(dir_plots, "/_raw/s6b_study_area.png"), width = 9, height = 18, dpi = 600, bg = "white")

p_s6b

Figure S7

Correlations of environmental variables

## extract covariates
## enough to use only 30 minutes data, because the others are subsets
df_boars_env <- 
  df_boars |> 
  filter(res == 30) |> 
  rename(t = "time") |> 
  with(track(x, y, t)) |> 
  extract_covariates(raster::stack(ras_env)) |> 
  mutate(hab = factor(hab, levels = c(0, 1), labels = c(" non-forest   ", " forest")))

df_boars_env$id <- df_boars |> filter(res == 30) |> mutate(id = paste("Wild boar", id)) |> pull(id)

## compute correlation matrix
cors <- 
  df_boars_env |> 
  mutate(hab = as.numeric(hab) - 1) |> 
  dplyr::select(
    "Forest\nhabitat" = hab, 
    "Imperviousness" = imp, 
    "Distance\nto water" = dtw
  ) |> 
  correlate(method = "spearman", diagonal = 1) |> 
  shave(upper = FALSE)

cors <- cors |>
 pivot_longer(
   cols = -term,
   names_to = "colname",
   values_to = "corr"
 ) |>
 mutate(
   rowname = fct_inorder(term),
   colname = fct_inorder(colname),
   lab = ifelse(!is.na(corr), sprintf("%1.3f", corr), "")
  )

ggplot(cors, aes(rowname, fct_rev(colname), fill = corr)) +
  geom_tile() +
  geom_text(aes(
    label = lab, color = abs(corr) < .75),
    family = "PT Mono", size = 7
  ) +
  coord_fixed(expand = FALSE) +
  scale_color_manual(values = c("white", "black"), guide = "none") +
  scale_fill_distiller(
    palette = "PuOr", na.value = "white",
    direction = 1, limits = c(-1, 1),
    name = "Spearman\nCorrelation"
  ) +
  guides(fill = guide_colorbar(barheight = unit(10, "lines"), barwidth = unit(.6, "lines"),
                               draw.ulim = FALSE, draw.llim = FALSE)) +
  labs(x = NULL, y = NULL) +
  theme(
    axis.text = element_text(family = "PT Sans Narrow", size = rel(1.2)),
    axis.text.y = element_text(hjust = 0),
    legend.position = "inside",
    legend.position.inside = c(.88, .75), 
    legend.title = element_text(margin = margin(1, 0, 15, 0)),
    legend.text = element_text(hjust = 1)
  )

ggsave(paste0(dir_plots, "/s7_corr_covar_emp.png"), width = 8, height = 7, dpi = 600, bg = "white")

Figure S11

Environmental variables per location

p_s16 <- 
  ggplot(df_boars_env, aes(dtw, imp, color = hab)) +
  geom_vline(xintercept = 0, color = "grey65", linetype = "23") +
  geom_hline(yintercept = 0, color = "grey65", linetype = "23") +
  geom_point(shape = 21, fill = "white", stroke = 1.2) +
  geom_point(color = "white") +
  geom_point(alpha = .1) + 
  facet_wrap(~id, ncol = 4) +
  scale_color_manual(values = c("#a3811e", "#228822"), name = "Forest habitat:") +
  labs(x = "Distance to water (scaled)", y = "Imperviousness (scaled)") +
  theme(plot.margin = margin(1, 10, 1, 1))

ggsave(paste0(dir_plots, "/s11_env_vars_move_by_id.png"), width = 15, height = 12, dpi = 600, bg = "white")

p_s16


Session Info
# git2r::repository() ## remove information to ensure double-blind review
sessionInfo()
R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Ventura 13.6.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] here_1.0.1          corrr_0.4.4         amt_0.3.0.0        
 [4] patchwork_1.2.0     colorspace_2.1-0    rcartocolor_2.1.1  
 [7] rnaturalearth_1.0.1 terra_1.7-78        sf_1.0-16          
[10] purrr_1.0.2         lubridate_1.9.3     ggtext_0.1.2       
[13] ggplot2_3.5.1       stringr_1.5.1       forcats_1.0.0      
[16] tidyr_1.3.1         dplyr_1.1.4         readr_2.1.5        

loaded via a namespace (and not attached):
 [1] Rdpack_2.6                    DBI_1.2.3                    
 [3] s2_1.1.6                      remotes_2.5.0                
 [5] rlang_1.1.4                   magrittr_2.0.3               
 [7] e1071_1.7-14                  compiler_4.4.2               
 [9] systemfonts_1.1.0             vctrs_0.6.5                  
[11] profvis_0.3.8                 wk_0.9.1                     
[13] pkgconfig_2.0.3               fastmap_1.1.1                
[15] backports_1.5.0               ellipsis_0.3.2               
[17] labeling_0.4.3                utf8_1.2.4                   
[19] promises_1.3.0                rmarkdown_2.26               
[21] markdown_1.12                 sessioninfo_1.2.2            
[23] tzdb_0.4.0                    ragg_1.3.2                   
[25] rnaturalearthhires_1.0.0.9000 xfun_0.49                    
[27] cachem_1.0.8                  jsonlite_1.8.8               
[29] later_1.3.2                   R6_2.5.1                     
[31] RColorBrewer_1.1-3            stringi_1.8.4                
[33] pkgload_1.3.4                 Rcpp_1.0.12                  
[35] knitr_1.46                    usethis_2.2.3                
[37] httpuv_1.6.15                 Matrix_1.7-1                 
[39] splines_4.4.2                 timechange_0.3.0             
[41] tidyselect_1.2.1              rstudioapi_0.16.0            
[43] yaml_2.3.8                    codetools_0.2-20             
[45] miniUI_0.1.1.1                pkgbuild_1.4.4               
[47] lattice_0.22-6                tibble_3.2.1                 
[49] shiny_1.8.1.1                 withr_3.0.0                  
[51] evaluate_0.23                 survival_3.7-0               
[53] units_0.8-5                   proxy_0.4-27                 
[55] urlchecker_1.0.1              xml2_1.3.6                   
[57] pillar_1.9.0                  KernSmooth_2.23-24           
[59] checkmate_2.3.1               generics_0.1.3               
[61] sp_2.1-4                      rprojroot_2.0.4              
[63] hms_1.1.3                     commonmark_1.9.1             
[65] munsell_0.5.1                 scales_1.3.0                 
[67] xtable_1.8-4                  class_7.3-22                 
[69] glue_1.8.0                    tools_4.4.2                  
[71] d6berlin_1.0.0                fs_1.6.4                     
[73] grid_4.4.2                    rbibutils_2.2.16             
[75] devtools_2.4.5                raster_3.6-26                
[77] cli_3.6.3                     textshaping_0.4.0            
[79] fansi_1.0.6                   d6_0.1.0.4                   
[81] gtable_0.3.5                  ggsci_3.1.0                  
[83] digest_0.6.35                 classInt_0.4-10              
[85] htmlwidgets_1.6.4             farver_2.1.1                 
[87] memoise_2.0.1                 htmltools_0.5.8.1            
[89] lifecycle_1.0.4               httr_1.4.7                   
[91] mime_0.12                     gridtext_0.1.5